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The two-dimensional S = 1/2 XY model is investigated with an extensive quantum Monte 
Carlo simulation. The helicity modulus is precisely estimated through a continuous-time loop 
algorithm for systems up to 128 x 128 near and below the critical temperature. The critical 
temperature is estimated as Tkt = 0.3427(2) J. The obtained estimates for the helicity modulus 
are well fitted by a scaling form derived from the Kosterlitz renormalization group equation. 
The validity of the Kosterlitz-Thouless theory for this model is confirmed. 
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§1. Introduction 

The XY model in two dimensions has been discussed 
in various contexts such as magnets with easy-plane 
anisotropy, superconductivity in a thin layer, granu- 
lar superconducting materials, and insulator-superfluid 
transition in He 4 systems. Naturally, a number of works 
were devoted to clarifying the nature of the phase tran- 
sition and the low-temperature phase of this model. 
Among them notable was a large-scale Monte Carlo sim- 
ulation by Ding and MakivicBB Based on computa- 
tion of the linear in-plane susceptibility and the corre- 
lation length at various temperatures, they concluded 
that a phase-transition takes place at Tkt = 0.350(4)ErEP 
or 0. 353(3), a 1 and that this transition is of Kosterlitz- 
Thouless (KT)q) type. However, it is technically dif- 
ficult to distinguish an exponential divergence from an 
algebraic one. Because of this difficulty, the validity of 
their conclusion, an the nature of the phase transition 
was questionedBtP For the same reason, in spite of their 
vetVpfixtensive Monte Carlo calculation, Gupta and Bail- 
lidJ'El 1 did not draw a definitive conclusion although their 
numerical results seemed to suggest that the phase tran- 
sition of the classical XY model is exactly what the KT 
theory predicts. 

It should be noted that the above-mentioned technical 
difficulty is partially due to the absence of finite-size- 
scaling type analysis which is a common and powerful 
tool in_most numerical studies. However, Solyom and 
ZimanEP pointed out that the straightforward applica- 
tion of the ordinary finite-size-scaling analysis leads to a 
wrong result not only when the power-law temperature 
dependence of the correlation length, which is wrong, 
is assumed but also when the correct exponential diver- 
gence is assumed. They studied the size dependence of 
the first excitation gap in the S = 1/2 anisotropic XX Z 
model in one dimension, which is exactly solvable and 
known to have a transition of the KT type at the anti- 
ferromagnetic isotropic point. They found that the exact 
estimates for the finite systems do not fit into the stan- 



dard form of the finite-size-scaling at the critical point. 

Therefore, to obtain a definitive answer to the question 
concerning the nature of the phase transition, we need 
to have a correct form for the system size dependence of 
quantities of interest. In previous studies, to our knowl- 
edge, a systematic study of such system size dependence 
has been missing for the quantum XY model. There are, 
however, some reports on classical XY models. Instead 
of using the_ordinary finite-size-scaling form, Weber and 
Minnhag^nEi used the Kosterlitz renormalization group 
equationt 1 ! for the data analysis in their study of the 
classical XY model in two dimensions. They verified 
the KT type phase transition by comparing the size de- 
pendence of the helicity modulus, T, at the critical tem- 
perature with the renormalization group flow along the 
critical line that converges to the KT fixed point with log- 
arithmically slow convergence. They observed not only 
that the computed helicity modulus exhibits the logarith- 
mic dependence on the length scale but also that even 
the pre-factor of this logarithmic term agreed with the 
predicted one. Following the same idea, OlssorO per- 
formed a more detailed analysis of the classical model 
with an extensive Monte Carlo simulation. He observed 
that the system-size dependence of the helicity modulus 
agreed with the form derived from the Kosterlitz renor- 
malization group equation below and above the critical 
temperature as well as right at the critical temperature. 



libit the univer- 
This quantity 



The helicity modulus is known to e 
sal jump at the critical temperature.! 
corresponds to the super-fluid density when the model 
is regarded as a Boson system with hard-£ores. In the 
world-line quantum Monte Carlo methodjif 1 the helic- 
ity modulus is represented as the fluctuation in the to- 
tal winding number of world-lines by the following equa- 
tionM 



T = (T/2)(W 2 ), 



(1) 



where W = (W x ,W y ) with W x (W y ) being the total 
winding number in the x (y) direction. It is difficult to 
measure this quantity by means of a conventional world- 
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line quantum Monte Carlo method because a conven- 
tional algorithm is not ergodic in that the winding num- 
ber is not allowed to vary. In principle, it is possible to 
make it ergodic by introducing additional global move- 
ments of world-lines. In practice, however, such addi- 
tional and in most cases "ad-hoc" global movements are 
seldom accepted and the resulting estimates of the he- 
licity modulus have large statistical errors. Therefore, 
MakivicuiP divided the whole system into a number of 
sub-systems and measured the winding number for each 
sub-system which may vary. It is not totally clear if 
this alternative way of measurement gives the same an- 
swer as the conventional one does. Another difficulty of 
the conventional Monte Carlo method is its long auto- 
correlation time near and below the critical temperature. 
These difficulties limited the accuracy and the precision 
of the Makivic's estimates of the helicity modulus and 
narrowed the accessible temperature range of simulation. 

In the present paper, we report some results of the 
quantum Monte Carlo simulation, pf -the 5=1/2 XY 
model using the loop algorithmic lialij) with both dis- 
crete and continuous imaginary time representations^ 
Similar results for smaller systems with discrete imag- 
inary ±Lme representation have been reported else- 
where.EH> The use of the loop algorithm eliminates both 
of the above-mentioned difficulties. There the number 
of particles as well as the winding number can Jan. 
At the same time, there are a number of reportaiS'tif 
on the drastic improvement in the auto-correlation time 
of the simulation by loop algorithms. We aim at a 
detailed and precise comparison between_the quantum 
XY model and the theory by KosterlitzEf through an 
accurate estimation of the thermal fluctuation in the 
winding number near and below the critical tempera- 
ture. We show that such an estimation allows us to 
examine a new scaling form different from the ordinary 
finite-size-scaling. In this scaling form, the distance from 
the critical point, i.e., K — -Kkt appears in the form of 
(K - K KT ){log{L/L )) 2 , in contrast to (K - K K t)L vt 
in the ordinary finite-size-scaling. At the same time, the 
quantity x = ((7r/2)W 2 ) — 2 scales as x\og(L/Lo) rather 
than x/L" with some exponent uj. This "scaling" form 
is consistent with Olsson's fitting functions. 

In § ^, we describe loop algorithms on discrete and 
continuous imaginary time. In § ^, the definition of he- 
licity modulus, its improved estimator and details of sim- 
ulations are described. In § ||, estimates of the helicity 
modulus are presented and we summarize the results in 
§| 

§2. Monte Carlo Algorithms 

2.1 World-line Monte Carlo method for quantum XY 
model 

The 5=1/2 quantum XY model is defined by the 
following Hamiltonian. 

6 = E % = - J £($$J + ( 2 ) 

(ij) (ij) 

where (ij) runs over all nearest-neighbor pairs on a 
square lattice. As for the spin operators, we use the 



convention in which (S?) 2 = 1/4 (p — x,y,z). We will 
take J as the unit of the energy scale in what follows. 

In world-line Monte Carlo methods for d dimensional 
quantum systems, we generate world-line configurations 
in the (d + 1) dimensional space-time with probability 
density proportional to the integrand of the Feynmann 
path integral. In a conventional quantum Monte Carlo 
simulation, it is customary to discretize the imaginary 
time in the path integral form to derive a model with 
classical degrees of freedom in (d + 1) dimensions. We 
then perform an ordinary Metropolis-type Monte Carlo 
simulation of this classical system. This standard pro- 
cedure is called Suzuki-Trotter (ST) decomposition. \3 
Due to the recent development of the continuous time 
representation,! 2 !? it is no longer necessary to discretize 
the imaginary time by the ST decomposition. Neverthe- 
less, we start with it because it is still the most natural 
way to explain the algorithm. In the actual computa- 
tion, we have used both of the discrete and continuous 
time algorithms. 

Our Hamiltonian can be considered as a sum of four 
sub-Hamiltonians: H = Ha + Hb + He + Hd- Each 
sub-Hamiltonian is a sum of operators commutable with 
each other, i.e., 

Hx= k ij (X = A,B,C,D) (3) 

(ij)ex 

where 

A = E odd column, j = i + e x }, (4) 

B = {{ij)\i E odd row, j = i + e y }, (5) 

C = E even column, j = i + e x }, (6) 

D = {{ij)\i E even row, j = i + e y }, (7) 

and e x (or e y ) is a unit lattice vector in the x (or y) 
direction. 

Using the ST decomposition, we transform the parti- 
tion function into the following form, 

Z*52I[w(S p ). (8) 

S V 

Here, p stands for a plaquette in the (d+ 1) dimensional 
space-time with two edges perpendicular and the other 
two parallel to the imaginary time axis. The space-time 
location of its left-bottom corner is given by (i,t) with 

{(4n)Ar if % E A 

(4n + l)Ar if ieB 

(4n + 2)Ar if i £ C ' 1 > 

(4n + 3)At if i E D 

where the imaginary time spacing, At = (3/m, is the unit 
of the discretization of the imaginary time. The number 
of steps m is called the Trotter number. The symbol S p 
in eq. (||) is the local spin configuration on a plaquette 
p and w(S p ) is the two-spin propagator defined below. 
The symbol S in eq. (^) stands for the spin configuration 
of the whole space-time or the union of all Sp s, i.e., 
S = {J p Sp. We denote the four states of spin pairs by 
1 =TT, 2 =U, 3 =|t and 4 =L|. Then, the two-spin 
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propagator can be written explicitly as 

w(S p ) 
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.(10) 



The symbols 5^ mltlal ) anc j S p nu,iL> stand for the local state 
of two corners at the bottom and the top edge of the 
plaquette p, respectively. The local state of the whole 
plaquette, S p , can be regarded as the combination of the 
two. In the 4x4 matrix representation of eq. (10), the 
column index corresponds to four possible initial state, 



-r(final) 



s, 



(initial) 
P 



, and the row index S p 



(final) 



2.2 Loop algorithm with discrete imaginary time 

In a loop algorithm, we first assign a graph G p to each 
plaquette. A graph specifies the set of local states from 
which new state S' p is chosen. In other words, G p imposes 
a restriction on the local state. We call it a "graph" be- 
cause such a restriction can be conveniently expressed 
pictorially by drawing a line connecting two sites whose 
relative orientation is not allowed to be changed under 
the restriction G p . A group of spins connected by a se- 
quence of such lines is called a cluster. In the present case 
each site is shared by two plaquettes and it is connected 
to another site in each of the two plaquettes. There- 
fore, a cluster is a loop in the present case. Hence the 
name of the algorithm. In general, however, we need to 
introduce branching or connecting more than two sites 
in a plaquette. For example, in the case of the Ising-like 
XX Z model, we need to include a graph in which all four 
corners are connected. Once we assign graphs to all pla- 
quettes, we change the spin variables by regarding each 
cluster as a block spin and flipping it with an appropri- 
ate probability. In what follows, we briefly describe how 
we should choose the probability for assigning a graph to 
each plaquette and the probability for changing the spin 
variables. 

First, the two-spin propagator or the local Boltzmann 
weight w(S p ) can be written as 



w ( S p) = 2^ w (Sp,G p ), 
w(S p ,G p ) = v(G p )A(S p ,G p ). 



(11) 



Here the function A(S P ,G P ) takes on the value 1 when 
S p is compatible with the graph G p , and it takes the 
value 0, otherwise. 

For the S = 1/2 quantum XY model, the graph 
weight, v(G p ), is defined in Table I. The summation 
in eq. (11) is taken over the three graphs depicted in 
Fig. 1. We represent A(S P ,G P ) by the matrices in the 
same manner as eq. (10): 

/ 



A(-,3) = 



/ 1 



V o 



\ 




1 / 



(12) 



A Markov process generated by a traditional world- 
line algorithm stays in the spin configuration space: 

5(1) ^ 5(2) _> 5 0) _... t (13) 

whereas in a Markov process of a loop algorithm the spin 
and the graph configuration space alternate: 

5(D ^ G (D ^ 5(2) G (2) ^ 5O) ^ G (3) _^ . . . 

(14) 

The transition probabilities between the spin and the 
graph configuration are defined as follows: 



P(S -» G) 



Y[ p w(S p ,G p ) 
E G < UpHSp,G' p y 



P(G 



g s = U p w{Sp,G p ) 



(15) 



(16) 



where G = (J G p is the global graph. It is easy to 
check that these transition probabilities satisfy the de- 
tailed balance. 

We see from eq. (16) that P(G — > S) is vanishing if 
S p is not compatible to G p for any p, while it takes on a 
value independent of S as long as all S p 's are compatible 
to Gp's. This means that sites belonging to the same 
loop should be flipped simultaneously with a probability 
1/2, and two distinct loops should be flipped indepen- 
dently. Since flipping a loop is a global update, one may 
imagine that the cluster algorithm reduces the autocorre- 
lation time due to slow relaxation modes associated with 
large structures. TJais-is indeed the case as reported in a 
number of articles.HOJj 

2.3 Loop algorithmwith continuous imaginary time 

Beard and WieseEJ pointed out that we can take the 
continuous imaginary time limit (At — » 0) of the above- 
mentioned algorithm to obtain another algorithm which 
works directly on continuous time, instead of comput- 
ing quantities using discrete time with various Ar's and 
extrapolating the results to get estimates in the contin- 
uous time limit. In the discrete time version, we need 
to be careful in decomposing the Boltzmann operator in 
order to keep the correction term of high order in At. 
This is why we have to split the Hamiltonian into several 
pieces as we did in previous sections (eq. @)- In deriv- 
ing the continuous time version, however, we can avoid 
this complication because any correction term vanishes 
in the continuous time limit. In other words, we have to 
keep only the lowest order term to get the correct algo- 
rithm. This fact makes the derivation of the algorithm a 
little simpler as we see below. 

The matrix elements of the partial Hamiltonian H^j 



A(-,l) = 
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1 / 
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(sf nal >|i^4 initial) > = - 



G v 



a(G p )A(Sp,G p ). (17) 



(For the 5=1/2 quantum XY model, the a(G p ) is 
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defined in Table I). Then, for small At, the weights 
w(S p ) are expressed, up to the lowest order, as 



w{S p )^I(S p ) + ATj2a{G p )A(S p ,G p ), 



(18) 



where I(S P ) is the matrix element of the identity opera- 
tor on a plaquette with a "height" At in the imaginary 
time direction. For simplicity, we set A(S P , 1) = I(S P ) in 
the following. Taking the limit of an infinitesimal imag- 
inary time spacing (At — > 0), transition probabilities 
(eq. (15)) on a plaquette are reduced to the followings: 
1. If the present state S p is compatible to the graph 
G p — 1, i.e., if there is no exchange of spins in the 
time interval (t, t + At) , 

P(S P - G p ) = (AT)a(G p )A(S p , G p ) (G p ± 1), 
P(S P - 1) = 1 - P ( S p - G 'p)- ( 19 ) 



2. If the present state S p is incompatible to the graph 
G p = 1, i.e., if the state at t is different from that 
at t + At, 



#1 °(Gp) ^p) 



P(S P 1) = 



(20) 



We apply eq. (19) to plaquettes (with infinitesi- 
mal "heights") at which world- lines continue along 
the Trotter direction. The probability of choosing a 
graph Gp ^ 1 in the imaginary time interval At is 
(AT)a(G p )A(S p ,Gp). This means that in the continu- 
ous limit we distribute graphs G p over such an interval 
uniformly, i.e., with a Poisson process, with the probabil- 
ity density a(G p )A(S p , G p ). On the other hand, eq. (20) 
gives the probability (not probability density) with which 
a graph is assigned to each point of time where the lo- 
cal state is changed or two neighboring world-lines are 
exchanged. 

Consequently, for general models, the loop algorithm 
with continuous imaginary time can be summarized as 
follows. For each pair of nearest neighbor world-lines, 

1. distribute graphs G p (^ 1) with a Poisson process 
with eq. (19) over every imaginary time interval in 
which world-lines are not interrupted, 

2. choose a graph G p with eq. (20) at each point of 
time where states are exchanged between the two 
sites, 

3. assign the graph G p — 1 elsewhere, 

and then update spin values by flipping clusters 
(cq. (16)). To restate this procedure for the specific case 
of S = 1/2 quantum XY model, 

1. for each uninterrupted time interval during which 
spins on these world-lines are antiparallel, gener- 
ate "horizontal reconnections" (graph G p — 2 in 
Fig. 1) of world-lines with probability density J/4, 
and for parallel spins, generate "diagonal reconnec- 
tions" (graph G p — 3) of world-lines with probabil- 
ity density J/4, 

2. at each point of time where states are exchanged, 
assign a horizontal or diagonal reconnection with 



equal probability (1/2), 
3. assign the graph G p — 1 elsewhere. 

§3. Details of the Computation 

3.1 Helicity modulus 

In the present paper, we focus on the helicity modu- 
lus because it may be directly related to the solution of 
the Kosterlitz renormalization group equation and there- 
fore may exhibit some behavior characteristic of the KT 
transition. We have measured the helicity modulus by 
computing the fluctuation in the total winding number 
of world- lines with eq. (|l]). The total winding number 
W x (or W y ) is defined as the sum of winding numbers of 
individual world-lines. The winding number of an indi- 
vidual world-line is the number of times the world-line 
wraps around the system in the x or y direction before 
coming back to its starting point. The total winding 
number of world-lines for up spins is exactly the same in 
magnitude as that for down spins and has the opposite 
sign. We here count only the winding numbers for up 
spins. Alternatively, W x can be defined as 



x p 

±a x (Sp)^ c i S h 



(21) 



*Gp 



where L x is the lattice size in the x direction. The symbol 
OL x {Sp) stands for the function which takes on the value 
1 (or —1) if an up-spin world-line passes through the 
plaquette p in the positive (or negative) x direction and 
takes on the value 0, otherwise. Equivalently, c| is a 
constant which only depends on site i and takes on the 
value ±1/(4L X ) or zero. The other winding number W y 
can be also calculated in the same manner. 

3. 2 Improved estimator of helicity modulus 

In cluster type algorithms, it is often advantageous to 
re-express physical quantities in terms of graph variables 
rather than spin variables. The best known example is 
the graphical estimator for the linear magnetic suscepti- 
bility in the Swendsen-Wang algorithm for the uniform 
Ising model where the susceptibility is measured as the 
average cluster size rather than as the second moment of 
the sum of all spin variables. For the winding number, 
too, we can use a similar improved estimator. In our 
simulation, we have used this improved estimator for the 
squared winding number of world-lines. Equation (21) 
defines the ordinary estimator in terms of spin variables 
and it can be rewritten as 



W 2 = Wl + 



w x 



Wy 



(22) 
(23) 



where c| and are constants which only depend on 



site i. To be more specific, <f- = 1/4L X when the site i 
locates at the lower-left or upper-right corner of a shaded 
plaquette and c| = — 1/4L X otherwise. (Which is 'left' is 
an irrelevant question here as it is a matter of the overall 
sign of the winding number and we are interested only 
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in its squared value) . 



Here, we should notice that the winding number W x 

S) 



(0 



or Wy of a loop I formed in the process of the graph 
assignment can be expressed in terms of c| and c| defined 
above as follows, 

W«=2^c|5|, W«=2^c^|. (24) 
iei 



iei 



Then, the ordinary estimator O for the squared winding 
number can be decomposed into two parts: 



= + O 
where Oi mpr and O rem are defined as 



O, 



(25) 



(26) 



1 



0rem = \ fewpwp • ( 27 ) 

Since loops are flipped independently, the expecta- 
tion values of all the cross terms in O rem are vanishing. 
Therefore, 



{(D rern ) 0. 



(28) 



For the same reason, we can derive another useful equa- 
tion, 



(O 



-impr 



l )=o. 



Equation ( |2^ ) leads to 



(29) 



(30) 



This implies that the Oi mpr is another estimator of the 
squared winding number. The new estimator Oi mpr de- 
pends only on the graph variables. We can see that this 
new estimator is really "improved" as follows. 

Intuitively, one graph configuration represents 2 Nl spin 
configurations, where the Ni is the number of loops in 
the system. Therefore, a sampling of Oi mpr corresponds 
to taking an averaged value over many samplings of O. 
From eq. (p9|), the variance of O is related to those of 

Oimpr and Orem aS 

Var (O) = ((Oimpr -\- O r em) ) (Oimpr ~t~ Orem) v^-l) 



= Var (O impr ) + Var {O rem ) 

> Var (Oim pr ) ■ 



(32) 
(33) 



Equation (33) is valid for other improved estimators of 
the observable such as the one for the susceptibility or 
structure factor. Although we have introduced the im- 
proved estimator based on the discrete time representa- 
tion, we can use the same definition in the continuous 
time representation as well. 

We compared the error between two estimators in long 
runs at various temperatures. We found that the new 
estimator reduces errors in about twenty per cent and 
the total performance of the new estimator is about 1.5 
times better than that of the conventional one in terms 
of the computational time required. 



3. 3 Simulations 

In what follows, we present numerical results both 
from our older set of simulations performed with a dis- 
crete time version of the code and from our newer set 
with a continuous time version. We have taken various 
temperatures between 0.22 and 0.60 and used lattices 
with L = 8,12,16,24,32,48,64,96 and 128 in our sim- 
ulation. For simulations with the discrete time version, 
we have used m = 8, 16 and 32 for the Trotter number. 
When the systematic error of Trotter discretization ex- 
ceeds the statistical error, we have reduced the system- 
atic error by the extrapolation to m — 00 using three 
different Trotter numbers. For L = 12,24,48,96 and 
128 at all temperatures and for all L's near the critical 
temperature, we have used the continuous time version. 

The length of a typical run on L = 128 at each tem- 
perature is more than 10 6 Monte Carlo sweeps (MCS). 
The most time-consuming part in the entire code is the 
cluster identification. It is a task of assigning each spin a 
number that specifies which cluster the spin belongs to. 
In doing this, we only use the information of the local 
connectivity. To make a good use of vector processors for 
this kind of task, we need to use a vectorizable algorithm. 
For the discrete version of the code, we adaated an effi- 
cient vectorized code following Mino's ideaO This idea 
is based on the "divide-and-conquer" strategy. In this 
strategy, we firstly divide the lattice into many small sub- 
lattices and identifies clusters in each sub-lattice neglect- 
ing the connectivity outside of the cluster. This process 
can be easily vectorized or parallelized because cluster 
identifications of different sub-lattices are independent of 
each other. We then combine two adjacent original sub- 
lattices to form a larger sub-lattice and identify clusters 
in this sub-lattice. In doing this, we can use the informa- 
tion of the clusters in the smaller sub-lattices which have 
already been obtained in the previous step. We repeat 
this procedure until all the sub-lattices are combined into 
a single lattice, i.e., the original whole lattice. Using this 
algorithm, we achieved the efficiency of 1.5 million site 
updates per second per one vector processor on Fujitsu 
VPP500. For the loop algorithm on continuous imagi- 
nary time, we used a parallel computer and took trivial 
parallel approach. Our code does about 29000/ (L 2 /?) 
sweeps per second on a node of Hitachi SR2201. 

In our simulations, each run is divided into several 
bins. The length of a bin is taken large enough so that bin 
averages may be statistically independent of each other 
at least approximately. In order to check this condition, 
we have measured autocorrelation times of the improved 
estimator of the squared winding number at low temper- 
atures T < 0.35J by the standard binning analysis for 
a run of 10 5 MCS with L = 64. They turned out to be 
smaller than 2 MCS in all cases. We have not observed 
any difficulty due to low temperature except that, triv- 
ially, the computational time per one Monte Carlo step 
increases proportional to the inverse temperature. The 
statistical independence among bins is assured, because 
the smallest bin length used is 1000 MCS. An error bar 
shown in the figures in the present paper represents one 
standard deviation. Results of the squared winding num- 
ber (W 2 ) are summarized in Table II. 
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§4. Universal Jump in the Helicity Modulus 

4-1 Kosterlitz renormalization group equations 

The helicity modulus can be regarded as the renor- 
malized coupling constant that appears in the Kosterlitz 
renormalizatiorL_group equations. Furthermore, Weber 
and MinnhagerJlj) regarded the renormalization group 
flow of the solution as the system size dependence of the 
helicity modulus. They observed that the estimated val- 
ues at the critical temperature agree well with the theo- 
retical prediction derived from this idea. In the present 
paper, as we see below, we follow their idea but use a 
different method for the analysis in which not only the 
data at the critical temperature but also off-critical data 
are taken into account simultaneously 

The Kosterlitz renormalization group equations are 

dx o dy , „ 

-di = -y > i=~ xy - ^ 

Here, x and y are renormalized parameters after a renor- 
malization operation up to the length scale L = Lo(T)e l 
where Lq(T) is some characteristic length of the order of 
the lattice constant and has no singularity at T = Tkt- 
The renormalized coupling constant x is related to the 
helicity modulus and-the squared winding number by the 
following equations:H 



T 2 X ' 

Equation (S3) has an integral 



2. 



A(T)=x 2 (l)-y 2 (l), 



(35) 



(36) 



which does not depend on I = \og(L / L (T)) . As a func- 
tion of T, this integral has no singularity at T — Tkt- 
Therefore we can expand it in terms of the distance from 
the critical temperature, i.e., A(K) = a(K — -Kkt) + 
b(K - K KT ) 2 + ■■■ where K = J/T. The solution of 
cq. ( |34| ) is given by 



x{T,L) = 




(K > K KT ) 
(K = K KT ) 
(K < Kkt). 



(37) 



This solution is a special case of the following form: 



x(T,L) = r 1 f(Al 2 ). 



(38) 



This "finite-size-scaling" form can be obtained from the 
ordinary one 



x(T,L) = L UJ g (AL 1 /" 



(39) 



by replacing L by I = log(L / L (T)) . However, it is 
easy to see that eq. ( |38| ) cannot be made consistent with 
eq. ( |39] ) no matter what values are used for the expo- 
nents w and v. It should be also pointed out that even 
if we consider a little more general form for the ordinary 
finite size scaling 



x(T,L)=L«g{li{T)/L) 



(40) 



allowing £(T) to have the correct temperature depen- 
dence, it is still impossible to cast the solution eq. ( |37| ) 
into this form. 



From the solution eq. (|37j), it can be seen that the 
"scaling function" f(X) has no singularity at X=0. The 
scaling function should take on the value 1 at the critical 
point: 

f(X) = l + 0(X). (41) 

Note also that our scaling form eq. (^) is consistent with 
Olsson's fitting function [Equations (lla-c) with (16) and 
(18) in ref.|l|. 

4-2 Finite- size- scaling around T c 

In Fig. 2, the squared winding number is plotted 
against the temperature. We can see the strong system 
size dependence characteristic to the KT transition es- 
pecially around and above the critical temperature. Be- 
cause of this large size dependence it is virtually impos- 
sible to estimate the critical temperature and the criti- 
cal indices without knowing the "scaling form" that de- 
scribes the system size dependence. Figure 3 shows the 
same data rescaled using eq. (|8|). The parameters K KT 
and Lq are chosen to minimize the cost function defined 
in Appendix. 

In Fig. 4, the contour plot of the cost function is shown. 
The cost function is the measure of the "badness" of the 
scaling plot. It is basically the deviation from the local 
linear approximation. In order to eliminate data out 
of critical region, we have to select data points for the 
analysis. We have eliminated data points outside of the 
region, 



!<(w*><i + i, 

7T 7T 2 



(42) 



-1.5 < x = (K - K KT ){log{L/L )) 2 < 1.5. (43) 

The value of the cost function at the optimal choice of the 
parameters tends to be smaller as we eliminate more data 
points away from the critical point, making the resulting 
estimate more reliable. We have also selected data points 
with respect to the system size. However, if we eliminate 
too many data points the cost function does not have a 
meaningful minimum. As long as we obtain a meaningful 
minimum, the results do not significantly depend on the 
minimum system size adopted, as can be seen in Fig. 4. 
For the upper contour plot in Fig. 4 we have used data 
for system size L > 8, whereas for the lower we used 
L > 12. 

We here quote the estimates adopted in Fig. 3: 



T KT = 0.3427(2) J, L = 0.29(2). 



(44) 



The minimum value of the cost function turns out to be 
1.9. The deviation of this value from unity suggests the 
presence of correction to scaling. We can see from Fig. 3, 
that the value of the scaling function f(x) at x = is 
close to unity in agreement with the prediction (eq. (|4l"|)). 
The error in f(0) was determined in a similar fashion to 
those for Tkt and Lq. It results in 



/(0) = 1.0(1). 



(45) 



This agreement can hardly be explained unless the 
Kosterlitz-Thouless picture is assumed and is a clear in- 
dication for its validity in the present model. 
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4-3 Finite-size-scaling at T c 

We have alsoJried another analysis following Weber 
and Minnhagen.EEP namely, measuring the chi-square val- 
ues of the fitting to the critical form (The second equa- 
tion in eq. (|37])) as a function of the temperature. To be 
more specific, we assumed at each temperature the fol- 
lowing system size dependence of the helicity modulus. 

g^(W^A(T)( 1+ 21os(L ; Lo(T)) ) (46) 

This fitting form is expected to be correct only at the 
critical point with A(T KT ) = 1 (eq. (||) and eq. @). 

Since the number of data at each temperature is not 
enough, the critical temperature cannot be determined 
as the one at which the fitting is the best. Instead, we 
have tried two other procedures. One is to fix A(T) to be 
1 (but keep Lq as a fitting variable) and measure the chi- 
square values of the fitting. The result is shown in Fig. 5. 
From this figure we conclude that Tkt = 0.3430(5). The 
other procedure is to allow the coefficient A{T) to vary 
and see where A{T) crosses the line of A = 1, which 
should be the case at the critical temperature. The es- 
timated coefficient A(T) as a function of temperature 
while both A and Lq are allowed to vary is shown in 
Fig. 6. The critical temperature is estimated by a linear 
fitting of the data, yielding 



T KT = 0.34271(5)7, 
which is consistent with eq. (fl4|). 



(47) 



§5. Conclusion 

We have obtained accurate estimates of the helicity 
modulus as a function of temperature and system size. 
The results fit into the special finite size scaling form 
derived from the Kosterlitz renormalization group equa- 
tion identifying the renormalization scale with the sys- 
tem size. By this scheme we have avoided a technically 
difficult comparison between an exponential divergence 
and an algebraic one. The coincidence of the estimated 
critical value of the scaling function with the predicted 
one confirms the KT nature of the phase transition. In 
the numerical simulation we have used loop algorithms in 
both discrete and continuous time representations. Both 
of them have turned out to be quite efficient and ad- 
vantageous especially in estimating the helicity modu- 
lus which is usually a conserved quantity in conventional 
Monte Carlo simulation. 

A part of the numerical computation was done us- 
ing Fujitsu VPP500 at ISSP, the University of Tokyo 
and Kyoto University Data Processing Center and Hi- 
tachi SR2201 at the computer center of the University 
of Tokyo. N.K.'s work is supported by Grant-in- Aid for 
scientific research (No. 09740320) from the Ministry of 
Education, Science and Culture of Japan. 

Appendix: Cost Function for Evaluating Finite- 
Size-Scaling Plots 

In this appendix, we present a new estimator for evalu- 
ating finite-size-scaling plots. The estimator is basically 
a deviation from a local linear approximation analogous 



to the previous estimator ES> The advantage of this type 
of approaches is that we do not need to assume any spe- 
cific functional form for the scaling function. These es- 
timators are functions of trial values for various param- 
eters such as critical indices and temperature; Kkt and 
Lq in the present case. Since the previous estimator was 
discontinuous as a function of these parameters, search- 
ing for its minimum — the optimal point — was a rather 
tricky task. One could hardly tell if a local minimum 
found was really a global minimum. The discontinuity 
came from the reordering of the data points as we change 
the trial values of the exponents and the critical temper- 
ature. As we show in the following, we have eliminated 
this discontinuity by treating data points for different 
system size separately. 

In what follows, a data point consists of a rescaled 
estimate y, a rescaled error d and a rescaled parame- 
ter x at which y and d are estimated. In the present 
case, y = {{n/2){W 2 } - 2) log(X/L ) and x = (K - 
i^KT)(log(i/io)) 2 - Therefore, x, y and d depend on 
trial values for various parameters in general. Here we 
denote one standard deviation by d. Also we denote 
the i-th data point for the system size L by Pl,% = 
{zL,i,yL,ii ^L,i)- These data points are ordered so that 
< In the previous estimator, we dealt with 

all the data points without regard to the system size. 
Therefore, the order of the data points may change as 
the trial values for the exponents vary. Here we have no 
such reordering due to the change in the trial values. 

For each system size L', we connect adjacent points 
Pl 1 ,i and Pl'^+i by a straight line. Then, for each com- 
bination of a data point Pj,^ an d a system size L', we 
choose j such that xl',j < XL.i < xl',j+i- Next, we con- 
sider a point at x — xl,% on the line that connects two 
points, Pu.j and Pfj+i. We denote this point by Pfv 
Since the value y^ i is expressed as a linear combination 
of y L >j and yvj+i, 



y L L , 



_ XL',j + l - XL,: 



PVL',j + qVL',j+i, (A-l) 
q = l-p, (A-2) 

XL'j+l - %L',j 

the error d£ it which comes from errors dv j and dL',j+i, 
can be estimated through the ordinary error propagation 
rule. To be specific, 



P 2 d 2 



L',j 



(A-3) 



For convenience, we define i = Pl,%- In an ideal case 
where the linear approximation is good and the statis- 
tical and systematic errors are negligible, all Pf'/s for 
various L' should coincide with each other. 

Here we assume that the linear approximation is good 
and the systematic errors are negligible while the sta- 
tistical errors are not necessarily negligible. In order to 
evaluate the coincidence, then, we consider the weighted 
sum of deviations of the points P^a i- e -i 



(A-4) 
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with 

where n is the number of distinct system sizes. As statis- 
tics tells, the expectation value of Sl,i is 1. When we 
define S as the sum of all Sl/s, it evaluates quality of 
the finite-size-scaling plot. Its expectation value equals 
the number of data points N. The acceptable range of 
parameters is determined by 

S < mm S + m, (A-6) 

where m is of order 0(1). 



Table I. The v{G p ) and a{G p ) for S = 1/2 quantum XY model. 





v(G p ) 


a(G p ) 


1 


i(e-^ J + l) 


-y 


2 


±(e^-l) 


y 


3 


i(l-e"^) 


y 
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Table II. The squared winding number (W 2 ). 



Table II. The squared winding number (W 2 ) (continued). 



T/J 

/ 


L = 8 


L = 16 


L = 32 


L = 64 


0.220 




2.429(3) 


2.433(6) 


2.41(2) 


0.240 




2.214(3) 


2.225(5) 


2.22(2) 


0.260 




2.040(3) 


2.046(5) 


2.05(2) 


0.280 




1.876(2) 


1.865(4) 


1.86(2) 


0.300 




1.728(2) 


1.722(4) 


1.73(2) 


0.320 




1.592(2) 




1.57(1) 


0.325 


1.582(1) 


1.557(2) 


1.540(3) 


1.538(5) 


0.330 


1.552(1) 


1.524(2) 


1.508(3) 


1.493(3) 


0.332 


1.542(1) 


1.510(2) 


1.492(3) 


1.478(3) 


0.334 


1.528(1) 


1.492(2) 


1.478(3) 


1.471(3) 


0.335 


1.521(1) 


1.489(2) 


1.471(2) 


1.454(3) 


0.336 


1.518(1) 


1.483(1) 


1.463(2) 


1.450(2) 


0.337 


1.508(1) 


1.475(1) 


1.455(2) 


1.446(2) 


0.338 


1.505(1) 


1.466(1) 


1.444(2) 


1.433(2) 


0.339 


1.500(1) 


1.4617(9) 


1.442(1) 


1.425(1) 


0.340 


1.491(1) 


1.4556(9) 


1.434(1) 


1.418(1) 


0.341 


1.486(1) 


1.4480(9) 


1.426(1) 


1.407(1) 


0.342 


1.480(1) 


1.4417(9) 


1.419(1) 


1.400(1) 


0.343 


1.4727(9) 


1.4374(9) 


1.409(1) 


1.392(1) 


0.344 


1.4682(9) 


1.4292(8) 


1.403(1) 


1.385(1) 


0.345 


1.4633(9) 


1.4214(9) 


1.396(1) 


1.377(1) 


0.346 


1.4571(9) 


1.4149(9) 


1.388(1) 


1.366(1) 


0.347 


1.451(1) 


1.408(1) 


1.384(2) 


1.359(2) 


0.348 


1.445(1) 


1.401(1) 


1.370(2) 


1.348(2) 


0.349 


1.439(1) 


1.396(1) 


1.365(2) 


1.342(2) 


0.350 


1.433(1) 


1.391(1) 


1.356(2) 


1.334(2) 


0.351 


1.426(1) 


1.379(2) 


1.348(2) 


1.319(2) 


0.352 


1.423(1) 


1.376(1) 


1.338(2) 


1.309(2) 


0.353 


1.417(1) 


1.369(1) 


1.337(2) 


1.308(2) 


0.354 


1.412(1) 


1.360(1) 


1.327(2) 


1.298(3) 


0.355 


1.403(1) 


1.354(2) 


1.322(2) 


1.286(3) 


0.356 


1.395(2) 


1.349(2) 


1.311(2) 


1.280(5) 


0.357 


1.393(2) 


1.336(2) 


1.302(2) 


1.270(5) 


0.358 


1.386(2) 


1.336(2) 


1.298(2) 


1.262(5) 


0.359 


1.380(2) 


1.330(2) 


1.284(2) 


1.255(5) 


0.360 


1.373(1) 


1.322(2) 


1.279(2) 


1.247(4) 


0.362 


1.368(2) 


1.305(3) 


1.260(2) 


1.21(1) 


0.364 


1.356(2) 


1.296(2) 


1.246(2) 


1.20(1) 


0.366 


1.340(2) 


1.280(2) 


1.227(2) 


1.18(2) 


0.368 


1.328(2) 


1.262(2) 


1.207(2) 


1.14(1) 


0.370 


1.318(2) 


1.249(2) 






0.375 


1.290(2) 


1.214(2) 






0.380 


1.259(1) 


1.180(2) 






0.390 


1.203(1) 


1.105(2) 






0.400 


1.146(1) 


1.024(2) 


0.872(2) 


0.65(1) 


0.420 


1.033(1) 


0.856(2) 


0.598(2) 


0.242(6) 


440 


u. y 1 1 {I j 


n fi7fif9"\ 

U.U ( Ul ^ ) 






0.460 


0.804(1) 


0.492(2) 


0.142(1) 


0.007(1) 


0.480 


0.686(1) 


0.338(2) 


0.0518(7) 


0.0008(3) 


0.500 


0.579(1) 


0.214(1) 


0.0171(4) 




0.520 


0.479(1) 








0.540 


0.388(1) 








0.550 




0.0564(6) 






0.560 


0.310(1) 








0.580 


0.2465(9) 








0.600 


0.1935(9) 


0.0132(3) 







T/J L = 12 


L = 24 


L = 48 


L = 96 


L = 128 


0.339 1.475(1) 


1.449(1) 


1.432(1) 


1.416(1) 


1.416(1) 


0.340 1.468(1) 


1.442(1) 


1.423(1) 


1.411(1) 


1.409(1) 


0.341 1.462(1) 


1.432(1) 


1.414(1) 


1.401(1) 


1.398(1) 


0.342 1.4557(9) 


1.426(1) 


1.409(1) 


1.392(1) 


1.389(1) 


0.343 1.449(1) 


1.417(1) 


1.397(1) 


1.382(1) 


1.379(1) 


0.344 1.443(1) 


1.412(1) 


1.392(1) 


1.375(1) 


1.370(1) 


0.345 1.4368(9) 


1.404(1) 


1.383(1) 


1.366(1) 


1.363(1) 


0.346 1.4295(9) 


1.398(1) 


1.374(1) 


1.357(1) 


1.350(1) 
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Fig. 1. Three types of graphs for 5 = 1/2 quantum XY model. 
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Fig. 2. Hclicity modulus (or super fluid density) X = (T/2)(W 2 ) Fi §- 4 ' Contour P lots of the cost function for evaluating finite- 

as a function of temperature. The universal jump is expected at size-scaling plots. The plotted value is defined as S/N in Ap- 
the point where T = 2T/ir. Error bars are drown but most of 
them arc so small that they cannot be recognized. 



pendix. The smallest system size used is L = 8 for the top figure 
and L = 12 for the bottom. 




Fig. 3. A rcscalcd plot of the winding number fluctuation. The 
inset is a enlarged view near the critical point. 



Fig. 5. Chi-square values of the fitting with A(T) fixed to be 1. 
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